library(logr)
log_open("log_A7.log")
log_code()
####################
#REPLICATION FILES: APPENDIX 7
#Article: "When Does Online Public Diplomacy Succeed? Evidence from China's ‘Wolf Warrior’ Diplomats"
#Authors: Daniel Mattingly and James Sundquist
#This Version: June 18, 2022

####################
# Description:
# This file recreates the analysis presented in Appendix 7.

# Structure of script:
# 1. Load all observations (main analysis excluded suspected duplicates)
# 2. Fit PCA on this full set of observations
# 3. Fit models (ITT estimates)
# 4. Adjust p-values for both main analysis and ITT analysis
# 5. Report p-values, adjusted p-values, and Bonferroni significance in a table

# Note: These procedures are applied only to the pre-registered hypotheses,
# not every statistical test reported in the article and appendices.

#The script was written with R version 4.1.2
#It was created and tested on Mac OS X (11.6.5).


library(tidyverse)
library(xtable)
library(rio)
setwd("~/Dropbox/Research/Diplomacy/Replication")

### 1. Loading data
allobs <- import("allobs.csv")

### 2. Fit PCA
# China government
dv.names = c("Q21", "Q22", "Q25", "Q28", "Q29")
factor.obj = princomp(~ Q21+Q22+Q25+Q28+Q29, data=allobs, cor=TRUE, na.action = na.omit)
allobs$china_gov_pca = NA
allobs$china_gov_pca[complete.cases(allobs[,dv.names])] = as.vector(factor.obj$scores[,1])
allobs$china_gov_pca <- scale(allobs$china_gov_pca)
# China people
dv.names = c("Q17", "Q18", "Q24")
factor.obj = princomp(~ Q17+Q18+Q24, data=allobs, cor=TRUE, na.action = na.omit)
allobs$china_people_pca = NA
allobs$china_people_pca[complete.cases(allobs[,dv.names])] = as.vector(factor.obj$scores[,1])
allobs$china_people_pca <- scale(allobs$china_people_pca)
# China policy
dv.names = c("Q19","Q20","Q26_1","Q26_2","Q26_3","Q30")
factor.obj = princomp(~Q19+Q20+Q26_1+Q26_2+Q26_3+Q30, data=allobs, cor=TRUE, na.action = na.omit)
allobs$china_policy_pca = NA
allobs$china_policy_pca[complete.cases(allobs[,dv.names])] = as.vector(factor.obj$scores[,1])
allobs$china_policy_pca <- scale(allobs$china_policy_pca)
# China COVID
dv.names = c("Q66", "Q67", "Q68")
factor.obj = princomp(~ Q66+Q67+Q68, data=allobs, cor=TRUE, na.action = na.omit)
allobs$china_covid_pca = NA
allobs$china_covid_pca[complete.cases(allobs[,dv.names])] = as.vector(factor.obj$scores[,1])
allobs$china_covid_pca <- scale(allobs$china_covid_pca)


allobs$Treatment <- factor(allobs$Treatment, levels = c("control", "prochina", "antius", "antichina"))

### 3. Fit models
m.gov.all <- lm(china_gov_pca ~ Treatment, data = allobs)
m.people.all <- lm(china_people_pca ~ Treatment, data = allobs)
m.policy.all <- lm(china_policy_pca ~ Treatment, data = allobs)
m.covid.all <- lm(china_covid_pca ~ Treatment, data = allobs)

# Extract p-values
summary(m.gov.all)
summary(m.people.all)
summary(m.policy.all)
summary(m.covid.all)

# P-values for main analysis come from tables
# (replication code in Replication_main_results.R)

### 4. Bonferroni correction of p-values
# Main analysis
p.adjust(c(.011, .952), "BH") # Gov
p.adjust(c(.038, .326), "BH") # People
p.adjust(c(.003, .743), "BH") # Policy
p.adjust(c(.0001, .769), "BH") # Covid

# ITT analysis
p.adjust(c(.013, .988), "BH") # Gov
p.adjust(c(.036, .409), "BH") # People
p.adjust(c(.021, .930), "BH") # Policy
p.adjust(c(.0001, .962), "BH") # Covid

### 5. Construct table

Ys <- c("Government", "Government", "People", "People", "Policy", "Policy", "COVID", "COVID")
Treats <- rep(c("Foreign Aid", "Wolf Warrior"), 4)
p <- c(.011, .952, .038, .326, .003, .743, "< .001", .769)
p.adj <- c(.022, .952, .075, .326, .005, .743, "< .001", .769)
p.itt <- c(.013, .988, .036, .409, .021, .930, "< .001", .962)
p.itt.adj <- c(.027, .988, .071, .409, .042, .930, "< .001", .962)

pvals <- cbind.data.frame(Ys, Treats, p, p.adj, p.itt, p.itt.adj)
colnames(pvals) <- c("Outcome", "Treatment", "p-value", "B-H correction", "ITT p-value", "ITT B-H correction")
pvals

a6 <- print(xtable(pvals))
write(a6, file = "TableA6.tex")

log_close()
